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Abstract 

In this paper we consider a model based on branching process theory for the proliferation and the dissemination 
network of T cells in the adaptive immune response. A multi-type Galton Watson branching process is assumed 
as the basic proliferation mechanism, associated to the migration of T cells of the different generations from 
the draining lymph node to the spleen and other lymphoid organs. Time recursion equations for the mean 
values and the covariance matrices of the the cell population counts are derived in all the compartments of the 
network model. Moreover, a normal approximation of the log-likelihood function of the cell relative frequencies 
is derived, which allows one to obtain estimates of both the probability parameters of the branching process 
and the migration rates in the various compartments of the network. 

1 Moment computation for the network compartments 

In this Section we derive the first and second order moments of cell counts in the different compartments of the 
network model. A proliferation mechanism following a memoryless branching process known in the literature as 
a multi-type Galton Watson process is adopted. Population mean values and covariance matrices are essential 
tools both for stochastic simulation of the network and for implementing the parameter estimation procedure 
described in Section S2. We start our development by considering the source node of our network, i.e., the 
draining lymph node. 
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Draining lymph node 

The classical multi-type Galton Watson branching process (MGW)is a prototypical branching process, rep¬ 
resenting the evolution of a population whose members reproduce and die subject to random laws (see the 
classical textbooks {1, 2)). The generation of a cell, which is defined as the number division steps the cells 
undergoes before its birth, represents the type of that cell. We will denote by At the sampling time (or time 
step) of the discrete time process describing the proliferation. The generic time point t = n At will be denoted 
by n. Let the state of the population in the draining lymph node at time n be represented by the vector 
T^drin) = [Zdrfi{n), Zdr,i{n), Zdr, 2 {n), ■ ■ ■, Zdr,p{n)], where Zdr,i{n) are discrete random variables given by the 
counts of cells after i divisions, being p the maximum generation considered. 

It is well known that a MGW process is a homogeneous vector Markov process Zdr{n), n S N, with the following 
properties (see {1, 2)): 

PI) each cell action is independent from others; 

P2) each cell offspring generates its own branching process; 

P3) at each time point, the cell has no memory about previous time steps. 

In the most simple setting, a cell of type i can make three probabilistic decisions during each time step At: 
DI) remain in the same generation i with probability Sp, 

D2) divide to generate two cells of type z -I- 1 with probability jp 
D3) die with probability a^. 

The transition from Zdr{n) to Zdr{n -I- 1) in one time step is regulated by the probability generating function 
(pgf): 


f(s) = [/o(s),/l(s),/ 2 (s),...,/p(s)]^ , (1) 

where 

S — [^0, -Si, S2, . . . , Sp] G 5 I I — 1 5 

and 

/i(s) = djSj-I-Oi, z = 0,1,... ,p - 1, (2) 

/p(s) = otp SpSp. (3) 

It is clear that the i — th component of f (s) expresses all the possible outcomes of a cell of type z in one time 
step. Letting 1 = this concept is expressed by the property that /i(l) = 1, which implies that 
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Q!i + + 7i = 1; * = 0,1, 2,... — 1 and, because we assume that cells of generation p can only die or survive, 

ap + Sp = 1. This assumption is due to the fact that in our experimental setup, the CFSE fluorescence intensity 
becomes negligible for generations higher than p. 


The moments of the MGW process can be expressed in terms of derivatives of fi{s) (see {1, 2)). Specifically, 
the mean value of type j cell count at the time point 1 starting from a single cell of type i at the time point 0, 
is given by 


Mdr,%j = E[Zdr,j{l) I Zdr(O) = 6^] = 


dMs) 


dsn 


i, j = 0, l,...,_p, 


( 4 ) 


S = 1 


where e^, 0 < f < p denotes a vector whose i + 1 — th component is 1 and whose other components are 0. Since 
our Markov process is stationary, we can easily derive the transition matrix G which maps cell 

counts at time n to cell counts at time n + 1: 


Mdr = 


^ So 270 0 

0 271 

0 0 (52 

: : 0 


0 0 0 


0 

0 

272 

0 

0 


0 


Sp —1 27 ^ 


0 

0 

0 

'p-1 


\ 


( 5 ) 


It is easy to check that the mean value of Zdr{n) conditional to a given state Zdr{n — 1) is: 


E[Zdr{n) I Zdr{ri - 1)] = Zdr{n - l)Mdr ■ 


( 6 ) 


The mean value pi dr is obtained by taking expectaction with respect to Zdr(n — 1) in ([5]): 


= [^J'dr,o{'^), ■ • . ,/rdr,p(«)] = E[Zdr{n)] = “ l)Mdr ■ 


( 7 ) 


For computing the covariance matrix of cell counts Zdriji), we start by calculating the conditional second order 
moment E[Z'^^{n)Zdr{n) \ Zdriji - 1)]: 


E[Z^^{n)Zdr{n) \ Zdr{n- 1 )] = 

E[Z'^^{n) I Zdr{n - l)]E[Zdr{n) \ Zdr{n - 1)] + Cov[Zdr{n) \ Zdr{n - 1)] = (8) 

- l)Zdrin - l)Mdr + ELo ^i^drAn - 1) , 

where V/ represents the one step covariance matrix for one cell present in the state Z(0) = e; (see (I))'- 


{yih 


_ dfijs) dfijs) ' 

dsidsj dsi dsj J 


(9) 
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(Vih = 


d'^Ms) ^ dfijs) _ dfi{s) 


d'^St dsi 


ds^ 


( 10 ) 


- S=1 


By taking expectation with respect to Zdr{n — 1) in ([S]) we get the second order moment S'^j.{n): 

p 

= E:[Z'^r{n)Zdr{n)] = - l)Mdr + ^ itJi^dr,i{n - 1) . 


( 11 ) 


1=0 


Then, we get the covariance matrix Sdrin) recursion by considering that ^dr{n) = S’^^{n) — /x£,(n)/x^^(n): 

p 

Sdr{n) = M^^Sdr{n - l)Mdr + ^ V;/Xdr.;(n - 1 ) ■ ( 12 ) 

1=0 

Now, we want to extend this basic model by allowing the presence of migration from the draining lymph 
node. To this purpose, we introduce the r.v.'s rji i = 1, 2,... ,p — 1 representing binomial random variables 
associated to the cells of different types. We assume that = 1 if the cell of type i migrates in a time step At 
and ryi = 0 if the cell of type i doesn’t migrate. We denote by and 1 — m-i, i=l,2,...,p—1 the probabilities 
of the two possible events, respectively. Notice that, we assume here that the naive T cells (type 0) and the 
highest generation cells(type p) have null migration probabilities, i.e., ruo = nip = 0. This is in accordance with 
biological knowledge on the process and the fact that the CFSE measuring equipment has limited resolution. 
Now, we introduce the pgf’s describing the MGW process with the presence of migration: 


/dr.i(s) = 5i(l - mi)si + 7 i(l - mi+i)sf+i + Oi, i = 0,1,... ,p - 1 

(13) 

/dr,p(s) = T Oip . 

Notice that the parameters a^, * = 0,1,... ,p — 1 in the first equation above have a meaning different from the 
‘mortality’ rate of the basic model (without migration). Here, they take into account the cumulative effect of 
death and migration of cells of the various generations. 

Hence, we can use (l4|) with the pgf’s in (1131) for computing the transition matrix of the process in our source 
node Mjir G Rp+I’P+i in the presence of migration: 


M, 


/ 6o 270 ( 1 -mi) 0 0 

0 5i{l — mi) 271 ( 1 —m 2 ) 0 

0 0 (52(1 —m2) 272(1 —m3) 


dr 


V 0 


0 0 

0 

(l-mp_i)^p_i 27 p_i 


(14) 
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Then, with a slight abuse of notation, we will replace the mean value and covariance recursive equations © 
and (na with the following expressions holding in the presence of migration: 




(15) 


P 

Sdr{n) = Mj^Sdrin - 1)1^^^ + ^ V/^dr,/(n - 1) , (16) 

1=0 

where the one time step covariance matrices Vi are again computed through ([5]) and (uni). 

Now, we can derive the conditional expected value of the migrating T cells from the draining lymph node 
by introducing the probability generating function for one cell of type I in one time step: 


fmigfiis) = 7omisf 

/mig,i(s) = SlmlSl+Jlml+lsf_^_J^ ,l = l,...,p-2 , 

l(s) — lUTp—iSp—i 

fmig,p — 0 

S = [so,Si,S 2 ,---,Sp] e C^'+Ms/I < 1 ,; = 0,1,... ,p . 
Then, we get the conditional expected value of migrating T cells at time point n as: 


(17) 


E[Zrmg{ri)\Zdr{n - 1)] = Zdrijl - 1)M 


mig 


(18) 


where again the transition matrix G is computed through ([5]) with the pgf’s in (fT7|l : 




^rmg 

( 0 27 omi 0 0 

0 rriiSi 271TO2 0 

0 0 171262 272TO3 

: 0 : ■■■ 

: : : 0 TOp_i<5p_i 0 

y 0 0 0 0 0 0 


0 


0 


(19) 


By taking expectation over Zdr{n — 1) in (IT^ . we get the mean value of migrating T cells G 

at time point n: 
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To compute the covariance matrix at time point n, Smigin) G r(p+1'P+i)^ we first derive the conditional and 
unconditional second order moments of the migrating cell counts: 


| Zdr(n 1)] — 

mLgH I Zdr(n - imz^^g(n) I Zar(n - 1)] + Cov[Z^,g(n) | Z^,(n - 1)] = (20) 

^mzg^dr(^ “ ^)^dr(n — l)Mmig + mig ,l Zdr ,l{n — 1) , 


where 'Vmig,i represents the one step covariance matrix of one cell of type I available for migration from the 
draining lymph node. These matrices can be computed again through equations (jH]) and (jlOp , with the pgf’s in 

(HZD. 

By taking expectation over Zdr{n — 1), we get: 


p 

^migi.’’^) ~ '^[^mig(^)^iTiig(^)] = ^mig®dr(^ ~ l)Mmig + ^ ', '^mig,lfJ'dr,l{n' ~ !)■ 

1=0 

Then, the covariance matrix is readily obtained: 


Smigin) = S*^^gin) - ti^,gin)ti^,gin) 

— l)l^m2g T ~^mig^llddr^lin 1) • 


( 21 ) 


Transfer compartment 

Assuming a time step equal to 4 hours and that T cells take approximately 12 hours to move from the draining 
lymph node to the distal compartments, the transfer compartment can be built by assembling 3 serial subcom¬ 
partments TRl, TR2 and TR3 (see Fig. [1]). Of course, any different choice would not impact on the structure 
of the model, which is very flexible in this respect. Clearly, a rough knowledge of the migration time and of the 
time step, i.e., the mean division time, is necessary for appropriately structuring the transfer compartment. 



Figure 1: Network model schematization. We assume our network model composed by draining lymph 
node(dr), transfer compartment with its subcompartments TRl, TR2 and TR3 where the T cell immigrate in 
the mesenteric (mes), iliac (il) distal lymph nodes and spleen (spl). 


We assume that migrating T cells Zmigin) at time point n are located in the TRl subcompartment. We 
denote by Zruiin) the T cells counts vector of the first subcompartment. Then, the mean value fJ^rmin) G 
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and covariance matrix STRi{n) € ]R(p+i>p+i) are: 


fJ-TRi{n) = E[ZTRi{n)] = 

(22) 

^TRl{n) = Smig{n) . 

(23) 


After one time step in the TRl compartment, the cells ZtihI'/t-) undergo a transition while moving to the TR2 
subcompartment, according to the pgf ftr(s) = [/tr, 0 ) ftr,i, ■ ■ ■, ftr,p], defined as: 


ftr,o{s) = 0 

/tr,i(s) = /*(s), i = l,2,...,p-1 

/tr-.p(s) = 0 . 


From these pgf’s we easily derive the transition matrix G 


/ 0 0 0 

0 271 


Mf, = 


0 0 (52 272 0 0 

: 0 : ■■■ ■■■ 0 

: : : 0 Sp-i 0 

0 0 0 0 0 0 

Hence the conditional mean value of cell counts in the TR2 compartment is given by: 


V 


/ 


(24) 


ZTR2{n) = E[ZTR2in)\ZTRi{n — 1)] = ZTRi{n — l)Mtr ■ (25) 

By taking expectation over ZrRiin — 1) in (1^ . we get: 


MTfl2(^) — E[ZTR2{n)] — — \)Mtr ■ (26) 

The covariance matrix STR 2 {n) G at time point n can be easily computed following the same steps 

as for Smig'- 


p 

^TR2in) = M.Jl.STRli'n — + '^^'VtrJfJ,TRl,l{n — 1) , (27) 

1=1 

where 'Vtr,i is obtained through © and (ITTIl) by using the pgf’s introduced in (1^ . 

Concerning the splitting of the cell population among the the distal lymph nodes and spleen, we model this 
phenomenon as contemporary to the transition to the third transfer subcompartment. Hence, at the the end of 


7 




the third transfer time step, we are able to compute the three distinct subpopulations in the subcompartment 
TR3 (see Fig. 2). This step is described in the next subsection together with the proliferation in the distal 
nodes. 

Distal lymph nodes and spleen 

The cell population evolution in the distal lymph nodes and spleen depends on the intrinsic proliferation pro¬ 
cess and the immigration flow from the transfer compartment. Since we assume that splitting of the mi¬ 
grating population takes place during the transition from TR2 to TR3, we can compute the three flows of 
cells directed to the distal compartments, directly at the level of the third transfer subcompartment. The 
flow from the TR2 compartment towards the distal nodes is divided into three distinct components named 
'^TR 3 ,spi{'n),'ZTR 3 ,ii{n),ZTR 3 ,mes(ji) to identify the number of cells addressed to the spleen, iliac and mesen¬ 
teric lymph nodes, respectively. We will denote by Pspi, Pmes, Pu, {pu + Pmes + Pspi = 1) the probabilities that a 
cell decides to move to the spleen, the mesenteric or the iliac node. At the same time, our model takes care of a 
proliferation transition during the ‘splitting’ time step, according to the MGW pgf modified to account for the 
splitting. In the remaining part of this section, we will make reference to the spleen, giving for granted that sim¬ 
ilar formulas hold for the other distal compartments. The transition function from ZTR 2 {n) to ZTR 3 ^spi{n + 1) 
is regulated by the following pgf: 


^TR3,spl{s) = [frRSspipis), fTR3spl,l{s), fTR3spl,2{s), . . . , fTR3spl,pis)f 


(28) 


where 

s ^2; • ■ • 5 ^ G , ^ 1 , 

and 

fTR3spl,o{s) = 0 

fTR3spPi{^^ — Pspl^i^i T Pspl^i^iJ^-i T (1 Pspl^i Pspl^i') ^ i ~ 1 

fTR3spl,p{^') — (1 Pspl^p) T Pspl^p^p • 

Then, the conditional mean value of ZTR 3 ,spiin) is given by: 

'^TR3,spl{'n) = E[LTR3,spl{n)\ZTR2{n — 1)] = ZTR 2 {n — l)MspKt,spi (30) 


where 


'^^split^spl — Pspl^^tr • 


( 31 ) 



Expectation over 7iTR2{n — 1) leads to: 


MTfl3,spi(^) — ~ l)]V[spHt,spi (32) 

By using the same machinery adopted before, we get the recursion for the cell counts covariance matrix: 

p 

STR3,spi{n) = M'^pm^^piSTR2{n - l)'Msput,spi + {split,spi),itJ‘TR2,i{n - 1) , (33) 

1=1 

where 'V(^sput,spi),i the one step covariance matrix for one cell of type Z, which can be derived by exploiting 

Now, we turn to the computation of the mean and the covariance matrix of cell counts in the spleen. Before 
doing this, we have to quantify the immigrating flow Xspi{n) in the spleen and its covariance matrix: 

Xspiin) = E[Zspi{n)\Z,pi{n - 1) = 0] = fJ^TR3,spiin) 

(34) 

Wsp/(n) = Cov[Zspi{n)\Zspi{n - 1) = 0] = STR 3 ,spi{n) . 

Now, we can compute the conditional expected value of cell counts in the spleen, by considering the MGW 
process with the immigration flow { 40 )'- 

’Zspiin) = E[Zspi{n)\Zspi{n - 1)] = Zspi{n - l)Mdist + Kpi{n) , (35) 

where Xldist = By taking expectation over Zspi{n — 1), we get: 

l^spiin) = fJ^spiin - + Kpi{n) . (36) 

Finally, we can compute the covariance matrix Sspi{n) through the usual machinery, by suitably taking into 
account the immigrating population adding to the standard MGW process: 

p 

Sspi(n) = Mj,gjSsp/(n - l)Mdist + ^ 'ytr,ifJ^spi,i{n - 1) + Wsp/(n) , (37) 

1=1 

where Wsp/(n) is given by (l34ll and (1^ . 

2 Normal approximation for the log likelihood function of the cell 
relative frequencies 

All in vivo GFSE experiments require to sacrifice animals to collect data. This means that measurements taken 
at different time points actually refer to different individuals. This fact introduces an inter-individual stochastic 
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variability in terms of T cell counts which need be addressed in our inference scheme. 

To take into account this issue, model identihcation has been performed by using relative frequencies as model 
variables instead of T cell counts. This choice mainly based on the results in {3 ), allows one to derive an 
explicit expression for the asymptotic log-likelihood function of the relative frequencies when dealing only with 
one draining lymph node(5-5). When dealing with the entire network, the inference problem becomes much 
more complicated. Also in our context, we take a relative frequency approach, which will allow us to derive a 
closed form expression for a normal approximation to the log-likelihood function. To this purpose, we introduce 
the vectors Z(n) = [Zdr(n), Z*p,(n), Z* (n), Z;;,g^(n)] G ]RP“>‘ and /x(n) = G 

where ptot = 4p -|-1 and Zsp/(»^), Z* (n), /x*;(n), fJ^tnesin) G represent the cell counts and the 

mean values vectors computed in the previous section, with the exclusion of the first component representing 
naive (type 0) T cells. Similarly, we introduce the covariance matrices and build up the 

overall cell counts covariance matrix as 


S(n) 


^ Sdrin) 0 0 0 

0 Stpiin) 0 0 

0 0 S:,(n) 0 




G iPtot) 


V 0 


0 0 ) 


(38) 


We highlight that the matrix S(n) is introduced to incorporate the measurements taken at the different lymph 
nodes and spleen in the likelihood function. The block diagonal form is justified by the fact that the measure¬ 
ments performed in the different model compartments are independent. In order to construct the likelihood 
function, let A(n) G be the vector random variable representing the relative frequencies at time point n. 
The i — th component of A(n) is defined as 


Ai(n) 


Zi{n) 

Uin)’ 


(39) 


where U{n) = Zi{n). Notice that X)r=i Ai(n) = 1. Let r(n) be a vector whose components are: 


ri{n) 


Pi{n) 

ErriV.(n) 


1 ) • ■ • : Ptot ■ 


(40) 


Following arguments similar to those used in (5), let us introduce the matrix A(n) whose entries are given by: 


(A(n))ii 

= ('S'(n))-/^(l-ri(n)), 

(41) 


= -(-5'(n))-/^(l-rj(n)). 

(42) 

hj 

= 1^2^ ^Ptot^ j- 
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Moreover, define 

D(n) = A(n)'^R(n)A(n), 

where R(n) is the correlation matrix of the cell counts, whose entries are given by 


{R{n))ij 


iiS{n))T{S{n))y/ 


Then, the mean and the covariance matrix of relative frequencies can be derived as: 


(43) 


(44) 


E[A{n)] = 

S(n) = D(n)/(ErriV.W)^ . 

Since we know that X)r=i ^*(’^) = the relative frequencies distribution is degenerate. To this purpose, 
considering any subset of relative frequencies in the ptot dimensional vector of relative frequencies and the 
related covariance matrix, named (with abuse of notation) A(n) = [Ai(n),..., A(pj_,j_i)(n)]’^ G and 

S(n) G we get 

E[A{n)] = [ri(n),...,r(p,„^_i)(n)]'^ 

(45) 

S(n) = D(n)/(E^rlV.(n))^ 

where D(n) is defined like in (14311 through suitably rearranged submatrices of A(n). The above equations allow 
us to construct the normal approximation to the log-likelihood function for our estimation problem. In fact, 
assume that c„ mice are sacrificed at the time point n. Denote by Cfc(’^)j k = 1,2,..., c^, the relative frequency 
measurements at the timepoint n, and by Sfe(n) the relative frequencies covariance matrix at the same time 
point. Then, the contribution to the normal approximation of the log likelihood function at time n is given by: 

n) = ln(27r) - ^ ^ log{det{i:k{n)) - ^ ^(C/c(^^) - r(’^))^(S/c(n)"^(Cfe(n) - r(n)) , (46) 

k^l k^l 

where r(n) and are functions of the probability parameter vector 


^ To ; 5 Tl 5 ■ • ■ ) ^p—1: Tp—1 : TTl^ , . . . , Tflp— \ , Pspl , Pil 5 Pmes\ 


iT 


Of course, if independent measurements are taken at different time points n = ni,... ,nz, then the global 
negative approximate log likelihood function is 

Z 

C{e) =-Y,L{9-n{) . (47) 

i=l 

Notice that actually the cost function defined in (H51) would represent the true log likelihood function if the 
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relative frequencies were normally distributed: this is the case (asymptotically) when dealing with the draining 
lymph node only (see {3-5)). In our case, the cost function (ITZl) represents a normal approximation of the neg¬ 
ative log likelihood function, because we have no guarantee on the asymptotic normality of relative frequencies. 
By minimizing the cost function function (I47|) . we get parameter estimates: 

0 = arg{min(£(e)} (48) 

s.t. 

0 < di < 1, 0 < 7i < 1, di -I- 7i < 1, z = 0,1,... ,p - 1 

0<TOi<l, i = 1,... ,p — 1 

0<Sp<l, 

0 ^ Pspl ^ 1 ; 0 Pil ^ 1 ; 0 ^ Pmes ^ 1 ; 

Pspl 4 ” Pil 4 ” Pmes — 1 • 
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